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ABSTRACT 

Simulated hydroclimate variability in millennium-length forced transient and control simulations from 
the ECHAM and the global Hamburg Ocean Primitive Equation (ECHO-G) coupled atmosphere-ocean 
general circulation model (AOGCM) is analyzed and compared to 1000 years of reconstructed Palmer 
drought severity index (PDSI) variability from the North American Drought Atlas (NAD A). The ability of 
the model to simulate megadroughts in the North American southwest is evaluated. (NASW: 25°-42.5°N, 
125°-105°W). Megadroughts in the ECHO-G AOGCM are found to be similar in duration and magnitude 
to those estimated from the NAD A. The droughts in the forced simulation are not, however, temporally 
synchronous with those in the paleoclimate record, nor are there significant differences between the 
drought features simulated in the forced and control runs. These results indicate that model-simulated 
megadroughts can result from internal variability of the modeled climate system rather than as a response to 
changes in exogenous forcings. Although the ECHO-G AOGCM is capable of simulating megadroughts 
through persistent La Nina-like conditions in the tropical Pacific, other mechanisms can produce similarly 
extreme NASW moisture anomalies in the model. In particular, the lack of low-frequency coherence be- 
tween NASW soil moisture and simulated modes of climate variability like the El Nino-Southern Oscil- 
lation, Pacific decadal oscillation, and Atlantic multidecadal oscillation during identified drought periods 
suggests that stochastic atmospheric variability can contribute significantly to the occurrence of simulated 
megadroughts in the NASW. These findings indicate that either an expanded paradigm is needed to un- 
derstand multidecadal hydroclimate variability in the NASW or AOGCMs may incorrectly simulate the 
strength and/or dynamics of the connection between NASW hydroclimate variability and the tropical 
Pacific. 
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1. Introduction 

A particularly stark feature of proxy-estimated mul- 
tidecadal hydroclimate variability in the North Ameri- 
can southwest (NASW: 25°-42.5°N, 125°-105°W) is the 
occurrence of so called megadroughts [for a review, see 
Cook et al. (2007)]. Although drought definitions vary, 
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a megadrought can be defined as a persistent period of 
drought conditions lasting decades to centuries. Proxy 
records indicate the presence of two century-scale mega- 
droughts during the last millennium in the Sierra Nevada 
(Stine 1994; Cook et al. 2010), as well as a series of mul- 
tidecadal droughts that impacted much of the NASW 
(Cook et al. 2007; Herweijer et al. 2007). Understanding 
the cause of these megadroughts is important because of 
the potential for similarly extreme drought periods to 
emerge in the future. The presence of such drought re- 
gimes in the past (Woodhouse and Overpeck 1998; Cook 
et al. 2010; Herweijer et al. 2007) is particularly sobering 
when considering the current vulnerability of the NASW 
water supply to hydroclimate change (e.g., Schlenker 
et al. 2007). 

Assessing the ability of atmosphere-ocean general 
circulation models (AOGCMs) to simulate multidecadal 
drought features, like megadroughts in the NASW, is 
critical because these same models are used to make 
twenty-first-century climate projections. For example, 
the ensemble of International Panel on Climate Change 
(IPCC) Fourth Assessment Report (AR4) models pro- 
ject widespread drying in the subtropics over the coming 
century (Pachauri and Reisinger 2007; Seager and Vecchi 
2010). While this has been established as a forced re- 
sponse to increasing greenhouse gas concentrations in 
the atmosphere, it is unclear how forced and internal 
variability contribute to persistent hydroclimate fea- 
tures like past megadroughts. Determining the relative 
contribution is necessary because future hydroclimate 
will be determined by both radiatively forced changes 
and interannual-to-multidecadal internal variability. Hy- 
droclimate projections and associated risk assessments, 
therefore, require that AOGCMs capture both forced 
change and the amplitude and character of internal var- 
iability. Furthermore, the ratio of internal to forced hy- 
droclimate variability in the NASW has consequences 
for the predictability of future hydroclimate, particu- 
larly hydroclimate change related to anthropogenic 
greenhouse gas forcing. 

Despite the importance of testing the validity of sim- 
ulated variability in NASW hydroclimate on decadal- 
to-centennial time scales, it is difficult to do so with the 
instrumental record alone, therefore necessitating the 
use of paleoclimate estimates of past variability as model 
targets. Paleo model-data comparisons are thus vital 
exercises for evaluation of future hydroclimate projec- 
tions. The approach and execution of such comparisons 
will be investigated in this paper using millennium- 
length simulations from the ECHAM and the global 
Hamburg Ocean Primitive Equation (ECHO-G) AOGCM 
(Gonzalez-Rouco et al. 2006) and the North American 
Drought Atlas (NADA) (Cook et al. 2007). 


A wealth of research has implicated tropical sea sur- 
face temperatures (SSTs) as the dominant forcing of 
drought in the NASW. Schubert et al. (2004b, a), for 
instance, simulated the 1930s Dust Bowl drought as a 
response to tropical Atlantic and Pacific SST anomalies. 
Seager et al. (2005) and Herweijer et al. (2006) sub- 
sequently reproduced all of the major droughts of the 
nineteenth and twentieth centuries using an atmospheric 
general circulation model (AGCM) forced with ob- 
served SSTs. Similarly, Seager et al. (2008a) simulated 
megadroughts during the Medieval Climate Anomaly 
(MCA) with an AGCM forced with SSTs estimated from 
a single tropical Pacific coral record (Cobb et al. 2003). 
These simulated megadroughts were analyzed by Burgman 
et al. (2010), who noted similarities between the global 
pattern of modeled MCA hydroclimate and the one 
estimated from paleoclimate proxies. Herweijer et al. 
(2007) further analyzed megadroughts in the paleoclimate 
record, employing tree-ring reconstructions from the 
NADA to compare modern droughts to the megadroughts 
of the MCA. They proposed that the well-documented 
El Nino-Southern Oscillation (ENSO)-NASW tele- 
connection of the modern period (e.g., Seager et al. 
2005; Herweijer et al. 2006) was the likely forcing of 
persistent drought during the MCA, with the difference 
in drought persistence arising from the duration of 
drought-favorable SST conditions in the tropical Pacific. 
Similar work from Graham et al. (2007), using multi- 
proxy and modeling methods, also implicates the tropi- 
cal Pacific along with Indian Ocean SSTs as the principal 
influences on MCA hydroclimate changes. More re- 
cently, Feng et al. (2008) and Oglesby et al. (2011) have 
suggested that the tropical Atlantic played a role in 
forcing the MCA megadroughts, while Cook et al. (2013) 
have argued for the importance of dust aerosol forcing 
on both the spatial character and persistence of droughts 
in North America during the 1930s Dust Bowl and MCA 
droughts. 

Research has also placed the tropical Pacific SST 
boundary forcing in the context of model projections 
of global warming. Cook et al. (2010), for instance, 
recognize that, although the IPCC AR4 models robustly 
predict a shift toward dry conditions in the NASW, there 
is no agreement on the future state of the tropical Pacific, 
despite the strong connection between ENSO and NASW 
hydroclimate. This is because predicted drying arises not 
by any change in the spatial patterns of tropical SSTs but 
rather by overall planetary warming (Seager et al. 2007). 
Seager and Vecchi (2010), however, note that models 
that simulate an increase in the tropical Pacific east-west 
SST gradient (i.e., a trend toward more La Nina-like 
conditions) produce more drying in the NASW than those 
that simulate a decrease in the gradient. Recent work 
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nevertheless has complicated the gradient picture by 
demonstrating a large degree of internal variability of the 
zonal gradient in AOGCMs at centennial time scales 
(Karnauskas et al. 2012). 

Despite the large collection of literature in related 
areas, there are few analyses of megadrought occurrences 
and characteristics in simulations using AOGCMs. 
Meehl and Hu (2006, hereafter MH06) used a 1000-yr 
control run from the National Center for Atmospheric 
Research (NCAR) Parallel Climate Model (PCM) 
fully coupled AOGCM and found drought features of 
comparable length to proxy-estimated megadroughts 
that are mechanistically linked to low-frequency vari- 
ability in tropical Pacific SSTs. Additionally, Hunt (2011) 
analyzes global multiyear drought and pluvial occur- 
rences in a 10 000-yr control run of the CSIRO AOGCM 
and finds that persistent hydroclimate features can result 
from internal climatic variability, with stochastic atmo- 
spheric variability playing an important role. 

The following study builds on the work of MH06, 
Hunt (2011), and Herweijer et al. (2007) but differs in 
that we analyze both a forced transient millennium- 
length simulation and a 1000-yr control run together 
with 1000 years of proxy-estimated drought conditions. 
Two principal questions are addressed: 1) is the model 
capable of producing megadroughts that are character- 
istic of the paleoclimate record and 2), if so, are these 
drought features the result of internal variability or do 
they have a forced component? Answering these ques- 
tions is fundamental to understanding megadrought 
dynamics and interpreting simulations of future hydro- 
climate variability, which are in turn essential for future 
water supply management, risk assessment and infra- 
structure development in the NASW. 

2. Data and methods 

a. Observed and paleoclimate data 

Reconstructed PDSI data are from the North Amer- 
ican Drought Atlas (NADA) version 2a: the full details 
of which can be found in Cook et al. (2007). The data are 
reconstructed on a 2.5° latitude X 2.5° longitude grid of 
summer [June-August (JJA)] average PDSI values for 
the United States, as well as southwestern Canada and 
northern Mexico (286 grid points in total). The summer 
PDSI is reconstructed from a network of 1854 annual 
tree-ring records using a nested point-by-point regression 
method to produce records of maximal length. Verifica- 
tion statistics indicate that all grid points for the chosen 
analysis period (1000-1989 CE) and region (25°-42.5°N, 
125°-105°W) are highly statistically significant (decadal 
values of multiple determination, reduction of error, and 


coefficient of efficiency are greater than 0.7 for a 33-yr 
verification interval over California and Nevada; Cook 
et al. 2010). We also use SST data from the Kaplan ex- 
tended SST V2 product, which is a 5° latitude X 5°- 
longitude gridded SST field for the period 1856-present 
(Kaplan et al. 1998). 

b. Model 

Model analyses are performed using output from the 
ECHO-G AOGCM that combines the ECHAM4 and 
Hamburg Ocean Primitive Equation global (HOPE-G) 
atmospheric and ocean models, respectively (Legutke 
and Voss 1999). The resolution of the atmosphere is 
T30 horizontal (3.75°) by 19 vertical levels, while the 
ocean resolution is 2.8° in the zonal direction (T42) with 
equatorial refinement in the meridional direction that 
varies from 2.8° latitude at the poles to 0.5° at the 
equator with 20 vertical levels. The model employs a 
time-invariant adjustment of heat and freshwater fluxes. 
We use model SSTs, 2-m surface air temperature (SAT), 
precipitation, evaporation, sea level pressure, and soil 
moisture. The ECHO-G soil moisture model compo- 
nent is a single-layer bucket model with reservoir ca- 
pacity varying based on soil type (Legutke and Voss 
1999). For our purposes herein, the SAT, precipitation, 
and soil moisture are regridded from their native reso- 
lution to an even 2.5° X 2.5° grid. 

We use two ECHO-G simulations. The first is a 
1000-yr control simulation (clipped to 989 yr to match 
the length of the forced simulation used herein) that is 
run with constant external forcing set to mid-twentieth- 
century conditions. The second simulation is the ERIK2 
forced transient run (Gonzalez-Rouco et al. 2006) span- 
ning 990 years (1000-1990 CE; note that our subsequent 
analyses are over 989 years, 1000-1989 CE, owing to 
the employed yearly averaging interval of October- 
September) and driven by an estimated suite of external 
forcing factors including radiative effects of volcanic 
aerosols, concentrations of atmospheric constituents, 
and solar irradiance (Zorita et al. 2005). The run was 
initialized with pre-twentieth-century conditions and 
spun up for 100 years to the historical forcing of 1000 
CE (Gonzalez-Rouco et al. 2006). 

Internal variability of 2-m SAT, sea level pressure 
(SLP), and precipitation in the ECHO-G control run 
was evaluated by Min ct al. (2005a), who demonstrate 
that the model is capable of producing overall observed 
variability for all three of these variables. In a compan- 
ion paper, Min et al. (2005b) addressed the model treat- 
ment of interannual to decadal-scale internal variability 
using the same control run. They found that ENSO in 
the ECHO-G model exhibits stronger than observed 
amplitude and is too frequent and regular, with an 
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excessive spectral peak at 2 yr and muted variability in 
the 3-9-yr range. Importantly, the model produces rea- 
sonable ENSO spatial structures and teleconnections 
(Min et al. 2005b). Collectively, the ECHO-G simulations 
have been extensively analyzed (e.g., von Storch et al. 
2004; Zorita et al. 2003; Stevens et al. 2007; Gonzalez- 
Rouco et al. 2009; Karnauskas et al. 2012). Thus, de- 
spite the limited resolution of the ECHO-G model, its 
well characterized and vetted performance, and the 
availability of both millennium-length forced and control 
runs, make the EHCO-G simulations a comprehensive 
and consistent framework for an initial assessment of 
coupled atmosphere-ocean dynamics on decadal time 
scales. 

c. Drought variables 

For historical reconstructed estimates of NASW hy- 
drological conditions, we use PDSI from the NADA. 
For modeled drought conditions we use soil moisture 
(normalized over the length of the simulation) from the 
ECHO-G forced and control model runs. We use an- 
nually averaged soil moisture from the model, while 
the paleo-estimated NADA PDSI represents a JJA 
average. There were three principal motivations for us- 
ing the yearly averaged soil moisture instead of the JJA 
average or the model derived PDSI. First, annual-mean 
soil moisture provides a wide enough temporal window 
for addressing megadroughts and their dynamical causes. 
In particular, droughts driven by winter precipitation in 
the ECHO-G model (as are common for the NASW) will 
not necessarily produce a strong hydroclimate signal in 
the summer (JJA) average soil moisture (the modeled 
soil moisture memory time scale is 4-5 months versus 
12-18 months for the PDSI). Second, soil moisture is the 
model variable of most direct physical relevance to 
drought. Third, the PDSI appears to have a potentially 
troublesome dependence on temperature that causes 
a strong drift toward negative values in the model during 
the twentieth century, which is in neither the NADA 
PDSI nor the model soil moisture. The negative drift of 
PDSI in the ECHO-G model is likely a consequence of the 
temperature dependence of the Thornthwaite equation 
for potential evaporation used in the PDSI calculation and 
the excessive contemporary temperature trends in the 
model. Similar issues with the temperature dependence of 
PDSI have been raised in the literature (Sheffield et al. 
2012; Burke and Brown 2008; Milly and Dunne 2010). 
Normalized soil moisture is thus chosen as it provides 
a comparable analog to the PDSI, which is intended to 
represent a locally normalized anomaly of moisture sup- 
ply and demand (Dai et al. 1998, 2004). 

Although the simple bucket soil moisture scheme in 
the ECHO-G model may be problematic, soil moisture 


captures the influence of temperature and precipita- 
tion variability (see the subsequent variable comparison 
section) and considers snow water storage and melt, 
water interception by vegetation during rain or snow- 
melt episodes (skin reservoir), and water infiltration and 
runoff (Diimenil and Todini 1992). Nevertheless, it must 
be noted that more complicated soil vegetation atmo- 
sphere transfer (SVAT) schemes have the potential to 
alter results. 

d. Drought indices 

A drought index was calculated for the NASW by 
spatially averaging the normalized grid point anomalies 
of soil moisture for the NASW region. This box is some- 
what more restricted than that of MH06 and Herweijer 
et al. (2007) in order to maintain a homogeneous sample 
area (as determined by analyses of the spatial variance 
of the soil moisture field in the forced and control model 
runs). 

Drought definitions vary in terms of both input data 
(e.g., PDSI versus precipitation in the observed record) 
and criteria. We employ a drought definition similar to 
that described in Herweijer et al. (2007), with a drought 
commencing after two consecutive years of negative soil 
moisture anomalies and continuing until two consecu- 
tive years of positive anomalies (2S2E). Herweijer et al. 
required one year to start a drought and included a cri- 
terion based on spatial extent, which is not used herein. 
The adopted definition is different but broadly consis- 
tent (see section 3) with the drought definition of MH06, 
who define drought as consistently negative anomalies 
in an 11-yr running mean time series of box average pre- 
cipitation (droughts begin with the first year of anomalously 
negative precipitation and end in the first year of anoma- 
lously positive precipitation in the filtered time series). 

Droughts identified using the 2S2E definition were 
ordered by creating a drought density rank. For each 
drought period, the NASW index was summed from the 
first to the last year of the drought. These values were 
subsequently ranked by the negative value of the sum. 
This drought density ranking was chosen over a purely 
length-based ranking in order to incorporate both the 
persistence and severity of each drought. 

e. Variable comparison 

The annually and spatially averaged normalized pre- 
cipitation, precipitation minus evaporation ( P — /:’), and 
soil moisture over the NASW region are highly corre- 
lated in the model simulations (e.g., there is a 0.86 cor- 
relation between the soil moisture and precipitation 
NASW indices). Furthermore, yearly averaged soil mois- 
ture closely resembles the JJA average soil moisture for 
the ECHO-G model with a correlation of 0.7 between 


1 October 2013 


COATS ET AL. 


7639 


(a) Forced Soil Moisture Box Avg. Anomaly for NASW 
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FIG. 1. NASW (25°M0°N, 125°-105°W) box averaged (a) unnormalized forced soil moisture 
(m) and (b) forced PDSI (nondimensional) for the period 1000-1989 from the ECHO-G 
simulations (standardized over the period 1000-1850). (c) The NADA PDSI (nondimensional) 
for the same NASW region (standardized over the period 1931-90 using instrumental PDSI). 
The soil moisture is a yearly average, while the PDSI is a JJA average (again to match the 
NADA). The five most severe droughts using the described classification and ranking methods 
in section 2c are highlighted in red, with the 20-yr low-pass filtered time series plotted in blue. 


the two indices for the control run. The use of yearly 
average soil moisture is further justified by the agree- 
ment between the droughts identified in the annual and 
JJA control soil moisture indices (8 of the 10 largest 
droughts in the control run are in agreement using the 
drought identification and ranking methodology out- 
lined in the above section for the annual and JJA soil 
moisture). 

We also calculated model PDSI to allow for a direct 
comparison between simulated soil moisture and PDSI 
variability in the ECHO-G model (following Cook et al. 
2013). Model PDSI is derived on an even 2.5° X 2.5° grid 
using simulated precipitation and surface temperature 
as inputs. At each grid point PDSI was calculated and 
then standardized against a preindustrial normalization 
period (1000-1850 CE). Soil moisture capacity was speci- 
fied as 25.4 and 127 mm in the top and bottom layers, re- 
spectively, and evapotranspiration was calculated using 
surface temperature via the Thornthwaite (1948) method. 
The evolution of JJA PDSI is found to be comparable 
to yearly average soil moisture with three out of the 
five identified most severe droughts in agreement for the 
forced run (as can be seen in Fig. 1). 

3. Results and discussion 

a. Analyses of drought indices 

Figure 1 compares the simulated soil moisture in the 
forced model run to the simulated PDSI and proxy 


reconstructed PDSI. The simulated NASW soil mois- 
ture variability in the ECHO-G model correlates well 
(r = 0.75) with calculated model PDSI during the 1000- 
1900 period. The two records diverge in the postindus- 
trial period owing to an unrealistically large negative 
twentieth-century PDSI trend in the model simulation 
that can be attributed to an excessively positive tem- 
perature trend — more than twice the observed trend — 
and slightly negative precipitation trend in the twentieth 
century. If the modern/postindustrial period is ne- 
glected, the identified droughts using the two variables 
are consistent in three out of five cases. The exceptions 
are the late 1500s drought, which is the least severe of 
the five droughts identified in the soil moisture time 
series, and the twelfth-century drought in the PDSI time 
series; this latter drought is present but much smaller in 
magnitude for the soil moisture time series (and thus not 
identified as one of the five largest droughts in the time 
series). The disagreements in drought timing and mag- 
nitude seem to be associated with strong temperature 
controls on calculated PDSI that are not reflected in the 
modeled soil moisture response. This can be observed 
most dramatically in the PDSI estimates for the twen- 
tieth century in the forced run. 

In terms of drought severity, the model exhibits ap- 
proximately as much interannual and longer time-scale 
PDSI variability in the NASW region as the proxy re- 
cord (see the bottom two panels of Fig. 1). Although the 
PDSI has been noted to be difficult to compare in an 
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Fig. 2. Normalized soil moisture anomalies (forced and control model runs) and PDSI 
(NADA) for the period 1000-1989 CE averaged over the NASW region (25°-42.5°N, 125°- 
105°W): (a) the control soil moisture index, (b) the forced soil moisture index, and (c) the 
NADA PDSI index. Annual anomalies (black lines) are shown along with smoothed versions 
using a 20-yr low-pass filter (blue lines). The red highlighted periods in the annual time series 
are the five largest droughts as determined by the 2S2E drought definition and the drought 
density ranking. The gray shaded regions are the five largest droughts determined by the MH06 
drought definition. Note that the drought in the sixth century of the control run is actually split 
into two droughts using the 2S2E drought definition. The gray shaded drought in the first 
century of the control run is thus the sixth largest drought using our drought definition, (d) The 
volcanic and solar forcing time series (W m~ 2 ) used in the forced ECHO-G run for comparison 
to forced and NADA drought timing. Note that the bottom three panels are for the 1000-1989 
CE period; the timing of the control run in (a) is arbitrary. 


absolute sense (Dai et al. 1998, 2004), the model mega- 
droughts appear similar in severity to those in the pa- 
leoclimate record. An analogous comparison between 
the forced and control simulations indicates that soil 
moisture variability is comparable in each. In particu- 
lar, both model simulations have the same soil moisture 
variance in the NASW (30.25 mm 2 ). Subsequent com- 
parisons of forced and control responses in normalized 
soil moisture time series, therefore, can be interpreted 
as equivalent in their range of variance and can be com- 
pared to the proxy-estimated PDSI time series. 

Figures 2a-c show the time series of the normalized 
soil moisture anomaly averaged in the NASW box for 
the forced and control simulations and the NADA PDSI. 
Using drought density, the five largest droughts were 
ranked by both the 2S2E (highlighted in red) and MH06 
negative running mean index (gray shaded regions). 
Note that the drought occurring near the end of the 


500th year of the control simulation is split into two 
distinct drought features using the 2S2E metric. Our 
definition, when ranked by drought density, is consistent 
with the MH06 definition for 12 out of the 15 droughts 
(4 out of 5 for each dataset). There are slight differences 
in defined length because filtering removes positive ex- 
cursions in the MH06 definition that would delay or end 
droughts using unfiltered data in our 2S2E drought 
definition. Despite this, it is clear that both drought 
definitions are identifying the largest negative excur- 
sions in the indices. Any discrepancies occur because of 
a ranking reversal of the fifth and sixth droughts (in the 
NADA and forced run) or the division of a persistent 
period of drought into two droughts (in the control run). 

b. Paleo model-data comparison 

There is little or no agreement in timing between 
droughts in the forced simulation and the NADA PDSI 
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Fig. 3. The black dots are the number of droughts in the NASW region in the forced (For) and control (Con) simulations and the NADA 
that are at least (from left to right) 10, 15, or 20 yr in duration. Box plots are determined from 1000 red-noise time series with the same 
statistics as the corresponding model or NADA indices (middle bar is mean, top and bottom bars are the 75th and 25th percentiles, and the 
whiskers are the full data range), (far right) The spectra using the multitaper method (Mann and Lees 1996) for the forced (red) and 
control (blue) soil moisture indices. The dashed lines are the 5th and 95th percentile confidence intervals for the forced multitaper 
spectrum. 


indices. There are, however, droughts in both control 
and forced runs that are characteristic of the proxy es- 
timates. In particular, the three time series in Figs. 2a-c 
demonstrate that megadroughts in both model runs are 
of comparable duration to those of the paleoclimate 
record. Although the model exhibits more positive ex- 
cursions during a given drought period in some cases, 
the average length of the five most severe forced and 
control-run droughts is approximately equal to that of 
the NADA estimates (19, 22, and 21 yr, respectively). 

The presence of droughts in the control run that are 
comparable in length and severity to the forced run 
suggests that internal climate variability can cause mega- 
droughts in the model. Although it is unclear if observed 
megadroughts are the result of radiative forcing, overlap 
between the forced model and proxy-estimated drought 
time series would be expected if the reconstructed forcing 
used to drive the model is realistic and the modeled 
megadroughts are a forced response. This is not the 
case. For instance, the low-pass correlation between 
the forced drought index and NADA PDSI index (0.023) 
is not significantly different from the range of low-pass 
correlations between the forced drought index and 1000 
red noise series with the same persistence as the NADA 
PDSI index (r = —0.014 and r = 0.075 are the 25th and 
75th percentiles, respectively). Furthermore, the control 
drought index is just as temporally synchronous with 
the NADA record as the forced drought index, also in- 
dicating that any overlap between the historical droughts 
and those in the forced run occur by chance. Finally, 
a direct comparison to the forcing time series can be 
made in Fig. 2d, which indicates that modeled mega- 
droughts do not have a preferred forcing state. For 


instance, the 1800s model drought occurs during a pe- 
riod of relatively low solar forcing and high volcanic 
activity while the 1300s and 1500s model droughts are 
contemporaneous with relatively high solar forcing and 
low volcanic activity. These results provide evidence 
that low-frequency NASW hydroclimate variability in 
the ECFIO-G simulations is not solely a response to 
radiative forcing changes. 

As a further line of inquiry, the numbers of droughts 
greater than a threshold length are plotted in Fig. 3. The 
model produces more droughts in each threshold length 
than the NADA record, but the number of droughts in 
the model and NADA fall within a narrow range. Also 
in Fig. 3, the droughts in each dataset are compared to 
those of 1000 red-noise time series with the same char- 
acteristics as the corresponding model or observation 
(i.e., first-order autoregressive coefficient, variance, and 
mean). Historical and modeled droughts are more pre- 
valent than those in the red-noise time series for longer 
time scales (greater than the 90th percentile for all three 
datasets for droughts of 15+ and 20+ yr) but not for the 
10-yr threshold (Fig. 3). This latter observation is not 
surprising as noise series with some persistence should 
be capable of producing periods of persistent negative 
anomalies. The greater drought persistence in the ob- 
served and modeled data for longer time scales nev- 
ertheless indicates that there are likely mechanisms 
creating persistence beyond first-order autoregressive 
variations that are responsible for megadrought occur- 
rences. Interestingly, the box plots indicate that there 
is more persistent drought in the control simulation than 
in the forced simulation. A comparison of the spectra of 
the control and forced drought indices (Fig. 3) suggests 
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Correlation Coefficient 


Fig. 4. Correlation coefficient maps between soil moisture (models) or PDSI (NAD A) NASW indices and SST fields, (from top to 
bottom) The correlation of the overlapping period of the NADA with the Kaplan SST dataset, the full control simulation, the forced 
simulation for the modern period (1857-1989 CE), and the forced simulation for the period 1000-1856 CE. (left) The full unprocessed 
data, (middle) the 10-yr low-pass-filtered data, and (right) the 10-yr high-pass filtered data. The forced run was split into two sections 
because a strong positive trend in eastern Pacific SSTs in the modern period (1870-1989 CE) coincides with a slightly negative trend in the 
forced soil moisture index and washes out the phase connection between the two fields. 


that the control run does in fact exhibit more power in 
the decadal to multidecadal range. 

c. Drought spatial patterns and teleconnections 

To investigate the influence of the tropical Pacific on 
drought variability in the NASW, we calculate the cor- 
relation of the yearly SST field with the NASW drought 
index: the former was averaged from May to April, and 
for the model output the latter was averaged from 
October to September to reflect a lag between the 
ENSO-driven precipitation anomaly and the soil mois- 
ture anomaly (the NADA PDSI is JJA average). These 


calculations were performed for the full period in the 
two model simulations and the 133-yr time overlap 
between the NADA and Kaplan SST datasets (1857- 
1989). Three analyses were completed: one with the raw 
data, one with the 10-yr low-passed data, and one with 
the high-pass filtered data (separated using a 10-point 
Butterworth filter). Results are shown in Fig. 4. The 
NASW region has a weaker connection to the tropical 
Pacific in the annual and high pass for both the forced 
and control runs than in the observational data (see 
Table 1 for the average correlation value in the Nino-3 
region). Despite the discrepancy, the model index is still 
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Table 1. Average correlation coefficients in the Nino-3 region 
between sea surface temperature and NASW box average soil 
moisture (ECHO-G) and PDSI (NADA). The full correlation- 
coefficient field is shown in Fig. 5. 



NAD A/Kaplan 
(1857-1989) 

Forced 

(1000-1856) 

Forced 

(1857-1989) 

Control 

Full 

0.420 

0.403 

0.449 

0.238 

High 

0.435 

0.444 

0.512 

0.266 

Low 

0.357 

0.155 

-0.030 

0.116 


highly correlated with the tropical Pacific. Furthermore, 
it captures the major spatial features of the observed 
correlation map, indicating that the model contains re- 
alistic though weaker teleconnections. 

The low-pass correlation map is relevant for the pur- 
pose of understanding what drives multidecadal drought 
variability. For the observations, the connection of NASW 
PDSI to the tropical Pacific is only slightly lower for low- 
frequency variations as compared to high-frequency 
variations. In the model simulations the control run 
maintains a connection to the tropical Pacific when low- 
pass filtered data are used (similar to MH06). The forced 


run, on the other hand, does not maintain this connec- 
tion; this results from a strong positive trend in eastern 
Pacific SSTs in the modern period (1870-1989) that co- 
incides with a slightly negative trend in the forced soil 
moisture index and washes out the phase connection 
between the two fields. With the modern period re- 
moved there is a moderately positive correlation for low 
frequencies in the tropical Pacific, but still much weaker 
than the observational record (the average correlation 
between NASW soil moisture and Nino-3 SSTs is 0.16 in 
the forced run versus 0.36 for the paleo-observed record). 
The frequency dependent relationships are further illus- 
trated in Fig. 5 in which the wavelet coherence of the 
NASW box average NADA PDSI and the Nino-3 index 
is shown for the full 133 yr of the instrumental period. 
Shown below the instrumental plot are wavelet co- 
herence spectra between three randomly selected 133-yr 
segments of soil moisture and the corresponding Nino- 
3 SST indices from the ECHO-G control run. As was 
seen in the correlation fields, the model clearly exhibits 
much less coherence in the decadal time range than the 
observations. Note that the low-pass filtered observa- 
tions also show a relationship between positive PDSI 


WTC: Kaplan Nino3-NADA NASW PDSI 



1876 1896 1916 1936 1956 1976 



Years 

Fig. 5. (top) Wavelet coherence of NASW box average PDSI from the NADA (JJA average) 
with Nino-3 box average SST (averaged May-Apr) over the common period 1857-1989 CE. 
(lower panels) Coherence of NASW average soil moisture (averaged Oct-Sep) and Nino-3 SST 
(averaged May-Apr) for three random 133-yr subsets of the control run. The arrows show the 
phasing direction, the colored contours show the magnitude of the coherence, and the black 
outline shows significance at the 95% level. Shaded regions are outside the significance win- 
dows. Note the higher coherence in the decadal range for the observed/proxy data in the top 
panel. 
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Fig. 6. Average DJF precipitation anomalies (mm day - *) for each of the five most extreme droughts ranked by drought density. Time- 
weighted composite averages for the forced and control simulations are also shown. Blue indicates above average precipitation and red 
below average precipitation. The square box is the NASW region. 


and cool Atlantic Ocean SSTs. Like the tropical Pacific 
correlation, this is much weaker in the model. 

d. Dynamical diagnostics 

Not surprisingly, given the climatology of the NASW, 
negative December-February average (DJF) precipita- 
tion anomalies are the dominant cause of the annual soil 
moisture signal during NASW droughts. Figure 6 shows 
maps of the DJF precipitation anomalies during each of 
the five megadroughts in the forced and control simu- 
lations, as well as composites over all droughts. The 
spatial features are consistent within each of the droughts 
and between the forced and control simulations, with 


a positive precipitation anomaly in the Northwest (for 
all but the 784-804 control drought) while the NASW is 
anomalously dry. This structure is reminiscent of a La 
Nina winter moisture anomaly resulting from a north- 
ward shift of the storm track (e.g., Sarachik and Cane 
2010 ). 

Figure 7 shows the forced and control equatorial Pa- 
cific zonal SST gradient index with the five largest 
drought periods identified in the corresponding NASW 
index highlighted in red [the gradient index was calcu- 
lated by taking the difference between SSTs averaged 
in a western equatorial box of 5°S-5°N, 150°E-160°W 
and an eastern equatorial box of 5°S-5°N, 130°-80°W 
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Years 


Fig. 7. DJF equatorial Pacific zonal SST gradient anomaly (°C) for the forced and control 
runs calculated using the zonal SST gradient index defined in Karnauskas et al. (2009). The 
five largest drought periods as determined from the NASW soil moisture index and by 
using the 2S2E and MH06 definitions are highlighted in red or gray shading, respectively. 
The sixth-century drought is actually two droughts that were split by the 2S2E drought 
definition. 


following Karnauskas et al. (2009)]. Considering the 
evidence for synchronous phasing between La Nina 
states and negative NASW soil moisture periods on both 
interannual and decadal-to-multidecadal time scales (in 
observations) one might expect the drought periods to 
be coincident with the largest positive excursions in the 
gradient index (the most La Nina-like state). This is not 
the case, however, and the state of the tropical Pacific 
does not appear to have a consistent and strong control 
over simulated low-frequency drought periods in the 
NASW (only the late-thirteenth- and twentieth-century 
forced droughts and the late-sixth-century control 
drought correspond to persistent La Nina states). 
Low-frequency ENSO variability is, therefore, not the 
only mechanism driving persistent moisture anoma- 
lies in the NASW in the ECHO-G model. Similar 
analyses of both the Pacific decadal oscillation (PDO), 
the leading principal component of monthly SST 
anomalies in the North Pacific Ocean poleward of 
10°N (Zhang et al. 1997), and the Atlantic multidecadal 
oscillation (AMO), the 10-yr running mean of Atlantic 
SST anomalies north of the equator (Enfield et al. 2001) 
indices, suggests that these oscillations exert a similarly 
weak influence on modeled NASW hydroclimate (both 
were analyzed as in Fig. 7; the low-pass correlation be- 
tween NASW drought and the PDO and AMO indices 
are given in Table 2). Furthermore, there is very little 
consistency outside of the NASW region in the seasonal 
and annual mean model fields of temperature and evap- 
oration during drought periods. By contrast, the winter 
half-year average [November-April (NDJFMA)] SLP 
field shows a high pressure anomaly over the North Pa- 
cific during nearly all of the megadroughts (Fig. 8). This is 
consistent with a northward shift of the storm track. For 


the forced simulation, the hemispherically symmetric 
SLP anomaly in the composite is reminiscent of La Nina, 
but the individual droughts tend not to exhibit charac- 
teristic ENSO-driven SLP symmetry. In the control run, 
the composite and individual drought patterns are even 
less characteristic of ENSO variability, suggesting that 
stochastic Northern Hemisphere atmospheric variabil- 
ity can drive persistent NASW drought in the model. 
Although the Arctic Oscillation (AO), defined as the 
leading mode of the monthly mean wintertime SLP 
following Thompson and Wallace (1998), is more tightly 
coupled to both the forced and control soil moisture 
index (Table 2) than the corresponding ENSO, PDO, or 
AMO indices, the correlation is still weak. It is therefore 
likely that the stochastic variability associated with the 
simulated decadal-length droughts in the NASW is 
the result of numerous atmospheric modes that are not 
readily described using a single atmospheric index. Given 
the very consistent spatial structure of the precipitation 
anomalies and the above characterization of SLP anom- 
alies, our collective analysis suggests that stochastic at- 
mospheric variability can produce persistent northward 
shifts of the storm track in the ECHO-G simulated cli- 
mate, similar to those seen during La Nina events, and 
thus drive megadrought occurrences in the model. 


Table 2. Ten-year low-pass correlation coefficients between the 
AMO. PDO, and AO indices and NASW box average soil moisture 
in the ECFIO-G forced and control simulations. 



AMO index 

PDO index 

AO index 

Forced 

-0.087 

0.052 

-0.212 

Control 

0.012 

-0.010 

-0.149 
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Fig. 8. Winter (Nov-Apr) SLP anomalies (mb) for each of five droughts identified using the 2S2E identification metric and the composites 

over all the drought years for the forced and control runs. 


4. Conclusions 

Megadroughts in the North American southwest, in 
forced and control simulations using the ECHO-G 
AOGCM, are similar in duration and magnitude to 
those seen in the paleoclimate record. The droughts in 
the forced simulation are not, however, temporally syn- 
chronous with those in the proxy record or the forcing 
time series, nor are there significant differences between 
the drought features simulated in the forced and the 
control runs. This indicates that model-simulated mega- 
droughts can result from internal variability of the 
modeled climate system, rather than as a response to 
changes in exogenous forcing variations. The frequency 
and persistence of megadroughts in the model and North 
American Drought Atlas suggests that mechanisms 


beyond first-order autoregressive variability are pro- 
ducing these drought features. Although the ECHO-G 
AOGCM is capable of simulating megadroughts through 
a persistent anomalous SST forcing in the tropical Pa- 
cific (e.g., one in the late-sixth-century drought in the 
control run and one in the late-thirteenth-century drought 
in the forced run), other mechanisms can produce sim- 
ilarly extreme moisture anomalies in the NASW in the 
model. In particular, the lack of low-frequency co- 
herence between NASW soil moisture and other mod- 
eled fields and the PDO, AMO, and AO indices during 
identified drought periods suggests that stochastic at- 
mospheric variability can contribute significantly to the 
occurrence of simulated megadroughts in the NASW. 
These results, while limited to a single model, demon- 
strate the importance of analyzing both forced and 
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control simulations in concert with the paleoclimate 
record. Stochastic variability has been shown to drive 
drought in models on interannual to decadal time scales, 
particularly in weakly teleconnected regions by Hunt 
(2011). In this instance, it seems plausible that stochastic 
atmospheric variability in the ECHO-G model can pro- 
duce storm track shifts (and associated hydroclimatic 
changes like NASW drought) that are uninterrupted by 
tropical Pacific influence because of the weak NASW- 
ENSO teleconnection on multidecadal time scales. 

In the observational record, persistent droughts in the 
NASW have all been tied to cool tropical Pacific SSTs 
(e.g., Seager et al. 2005; Herweijer et al. 2006), but it is 
not known if this relation holds for the entire last mil- 
lennium. Consequently, these model results have two 
implications depending on whether the modeled hydro- 
climate variability is a reasonable representation of the 
actual climate system: 1) if the model is accurately 
simulating real-world variability, then stochastic atmo- 
spheric variability and ENSO both appear capable of 
producing persistent droughts in the NASW or 2), if the 
model is misrepresenting the actual variability, then this 
feature is a likely component of AOGCMs that will in- 
fluence future projections of hydroclimate, an in- 
accuracy that must be addressed when assessing model 
projections. One possible explanation for point two is 
that a weak teleconnection between the NASW and 
the tropical Pacific Ocean in the model allows atmo- 
spheric variability to drive droughts, whereas the 
tighter link to the Pacific in nature ensures that mega- 
droughts are more strongly forced by tropical Pacific 
SST anomalies. Additionally, there is observational ev- 
idence that warm tropical Atlantic SSTs can create 
a tendency toward dry conditions in the NASW (Seager 
et al. 2008b; Kushnir et al. 2010; Nigam et al. 2011) and 
this has been appealed to as a cause of medieval mega- 
droughts (Feng et al. 2008; Oglesby et al. 2011). The 
connection of the NASW drought index in the model to 
the Atlantic is weaker than observed and this too could 
allow atmospheric variability to exert a stronger rel- 
ative influence on NASW hydroclimate. 

Longer records of proxy-estimated tropical Pacific 
SST (e.g., Emile-Geay et al. 2013) are necessary to 
assess the state of ENSO during megadroughts and to 
determine how coherent previous NASW drought and 
ENSO variability may have been prior to the observa- 
tional record. In the meantime, additional analyses of 
AOGCM simulations will identify what produces model- 
simulated megadroughts and help evaluate model treat- 
ment of regional low-frequency hydroclimate variability. 
In particular, a model intercomparison employing mul- 
tiple AOGCMs is necessary to determine if stochastic 
atmospheric variability similarly influences NASW 


megadrought occurrences in the most recent generation 
of AOGCMs. This will be possible as an increasing 
number of last-millennium simulations from the Paleo- 
climate Modelling Intercomparison Project 3 (PMIP3)/ 
phase 5 of the Coupled Model Intercomparison Project 
(CMIP5) archive (Taylor et al. 2011) become available. 
The analyses in this paper are thus a proposed frame- 
work for assessing model performance with regard 
to hydroclimate variability and for performing com- 
parisons to paleoclimate data in future multimodel 
comparisons. 
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